package gbench.sandbox.plot;

import org.junit.jupiter.api.Test;

import java.text.NumberFormat;
import java.util.ArrayList;
import java.util.stream.Collectors;

import gbench.common.fs.XlsFile;
import gbench.common.matlib.MatlibCanvas;
import gbench.common.matlib.MatlibCanvas.NVec;
import gbench.common.matlib.MatlibCanvas.Range;
import gbench.common.matlib.MatlibCanvas.Seq;
import gbench.common.matlib.MatlibCanvas.XVec;
import gbench.common.matlib.data.DataReader;
import gbench.common.matlib.img.ImageOps;
import gbench.common.matlib.img.ImageWrapper;

import static gbench.common.matlib.MatlibCanvas.Range.rG;
import static gbench.common.matlib.MatlibCanvas.println;
import static gbench.common.matlib.data.Pipe.DUMMY;
import static gbench.common.tree.LittleTree.Elem.cph2;
import static gbench.common.tree.LittleTree.IRecord.slidingS;
import static gbench.common.tree.LittleTree.LIST;
import static gbench.common.tree.LittleTree.ZIP;
import static java.lang.Math.PI;
import static java.lang.Math.random;

/**
 * 绘制多变量线性回归的损失函数的演变趋势
 * @author gbench
 *
 */
public class JunitPlotRegression {
    
    @Test
    public void foo() {
        final var nf = NumberFormat.getNumberInstance(); // 数值格式化器

        final var volume = 500; // 数据规模
        final var num_variables = 10; // 变量个数
        final var data = NVec.of(volume, i -> random()); //原始数据
        final var mx_pipe = DUMMY(Object[].class)
                .cast(Double[][].class)
                .bind(XVec::OF)
                .bind(xv -> xv.fmap(NVec::new)); // 把原始数据改造成 二维数组
        final var mm = data.reshape2(volume / num_variables, num_variables); // 二维数组矩阵
        final var x_train = mx_pipe.evaluate(mm); // 训练数据 X
        final var x_train_t = XVec.OF(Seq.transpose(mm)).fmap(NVec::new); // 训练数据 X 的转置
        println("x_train （训练数据 X）:\n", XlsFile.DataMatrix.simpleOf(x_train.fmap(e -> e.toArray()).toArray()));
        println("x_train_t （训练数据 X 的转置）:\n", XlsFile.DataMatrix.simpleOf(x_train_t.fmap(e -> e.toArray()).toArray()));

        final var actual_thetas = NVec.of(10, i -> i + 1);// 实际参数
        final var y_train = x_train.fmap(actual_thetas::dot).mutate(NVec::of); // 预期结果
        println("y_train（预期结果）:\n", y_train.toString(e -> e + "\n"));

        final var thetas = NVec.of(10, i -> 0); // 估计参数
        println("估计参数(初始状态,训练前:)",thetas);

        final var error_pipe = DUMMY((XVec<NVec>) null)
                .bind(xi -> xi.fmap(thetas::dot).mutate(NVec::new))// x_i_j*theta_j -> y_i
                .bind(y_i -> y_i.sub(y_train));// 偏差计算管道

        final var error = error_pipe.evaluate(x_train); // 计算偏差
        final var rate = 0.001; // 学习速率
        final var eps = 1e-5d; // 误差容忍度
        final var losses = new ArrayList<Double>(); // 损失函数

        var loss = 0d; // 偏差函数
        while ((loss = error.fmapN(e -> e * e).mean()) >= eps) { // 持续用估计参数计算 损失程度，直到达到容忍范围。
            thetas.update((j, theta_j) -> theta_j - rate * 2 * error.dot(x_train_t.get(j)));// 重新估计 thetas
            error.update(error_pipe.evaluate(x_train)); // 更新偏差计算结果
            losses.add(loss); // 计算本次的 损失程度。
        }// white

        // 输出估计参数
        println("估计参数(训练后)" + thetas.toString(nf::format));

        // 绘制损失的演变过程
        plot(losses);

        try { // 数据图的展示时间
            Thread.sleep(100000);
        } catch (InterruptedException e) {
            e.printStackTrace();
        }
    }

    /**
     * 简单的绘图
     *
     * @param points
     */
    public void plot(ArrayList<Double> points) {

        // 业务数据准备：x,y作为后缀 表示 业务（原始)数据
        final var dtx = NVec.of(points.size(), i -> i); // 业务数据空间的 水平采样点 DaTa X （补充水平采样点)
        final var dty = NVec.of(points); // 业务数据空间的 垂直采样点 DaTa Y

        // 绘图数据：x,y作为前缀的 表示 图像（绘图）数据
        final var xrg = Range.of2(100, 450);// x 的绘图区域范围 X RnGe
        final var yrg = Range.of2(50, 450);// y 的绘图区域范围 Y RnGe
        final var xtf = DataReader.transformer(xrg, dtx.rangeL()); // 水平坐标转换器 X TransFormer
        final var ytf = DataReader.transformer(yrg, dty.rangeL()); // 垂直坐标转换器 Y TransFormer
        final var xs = NVec.of(dtx.stream(i -> i * 70).map(xtf)); // 转换成绘图水平坐标的点 :X valueS，放大X100倍离散化数据， 以方便观察细节。
        final var ys = NVec.of(dty.stream().map(ytf)); // 转换成绘图垂直坐标的点 :Y valueS

        // 绘图操作
        final var size = 1000; //  绘图画布的x水平方向的数据范围尺寸
        final var matlib = new MatlibCanvas(size); // 创建数学家画布

        // 图形参数准备
        final var rhotetas = Range.cph(rG(0, 5).limit(10), rG(-PI, PI).limit(40))
                .collect(Collectors.toList()); // 小圆点的绘制向量长度rho与角度theta
        final var centers = LIST(ZIP("x,y", xs.list(), ys.list()).stream() // 组合成点结构的 圆心坐标序列
                .map(rec -> NVec.of(rec.toArray(Double.class))).map(matlib::xY));// cartesian 形式的 点 的 圆心位置 的 坐标 序列

        final var bimg = matlib.plot(size, size / 2, pen -> {
            slidingS(centers, 2, 1, true).forEach(pen::drawStraightLine);// 将序列剪裁收尾相连的选段序列，并连接成线。
            pen.drawPoints(cph2(rhotetas, centers).map(matlib::rhotheta_affine_offset)); // 绘制图像的点
        }, null, true, img -> new ImageWrapper(img, true).drawCoordinateSys(xrg, yrg, dtx, dty, xs, ys)); // 绘制并生成图片

        ImageOps.jpgwrite(bimg,"C:/Users/gbench/Pictures/绘制多变量线性回归的损失函数的演变趋势.jpg"); // 保存图片
    }

}
